home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / cblas / source_her2k.h < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  9.4 KB  |  270 lines
  1. /* blas/source_her2k_c.h
  2.  * 
  3.  * Copyright (C) 2001 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. {
  21.   INDEX i, j, k;
  22.   int uplo, trans;
  23.  
  24.   const BASE alpha_real = CONST_REAL0(alpha);
  25.   BASE alpha_imag = CONST_IMAG0(alpha);
  26.  
  27.   if (beta == 1.0 && ((alpha_real == 0.0 && alpha_imag == 0.0) || K == 0))
  28.     return;
  29.  
  30.   if (Order == CblasRowMajor) {
  31.     uplo = Uplo;
  32.     trans = Trans;
  33.   } else {
  34.     uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
  35.     trans = (Trans == CblasNoTrans) ? CblasConjTrans : CblasNoTrans;
  36.     alpha_imag *= -1;        /* conjugate alpha */
  37.   }
  38.  
  39.   /* form  C := beta*C */
  40.  
  41.   if (beta == 0.0) {
  42.     if (uplo == CblasUpper) {
  43.       for (i = 0; i < N; i++) {
  44.     for (j = i; j < N; j++) {
  45.       REAL(C, ldc * i + j) = 0.0;
  46.       IMAG(C, ldc * i + j) = 0.0;
  47.     }
  48.       }
  49.     } else {
  50.       for (i = 0; i < N; i++) {
  51.     for (j = 0; j <= i; j++) {
  52.       REAL(C, ldc * i + j) = 0.0;
  53.       IMAG(C, ldc * i + j) = 0.0;
  54.     }
  55.       }
  56.     }
  57.   } else if (beta != 1.0) {
  58.     if (uplo == CblasUpper) {
  59.       for (i = 0; i < N; i++) {
  60.     REAL(C, ldc * i + i) *= beta;
  61.     IMAG(C, ldc * i + i) = 0.0;
  62.     for (j = i + 1; j < N; j++) {
  63.       REAL(C, ldc * i + j) *= beta;
  64.       IMAG(C, ldc * i + j) *= beta;
  65.     }
  66.       }
  67.     } else {
  68.       for (i = 0; i < N; i++) {
  69.     for (j = 0; j < i; j++) {
  70.       REAL(C, ldc * i + j) *= beta;
  71.       IMAG(C, ldc * i + j) *= beta;
  72.     }
  73.     REAL(C, ldc * i + i) *= beta;
  74.     IMAG(C, ldc * i + i) = 0.0;
  75.       }
  76.     }
  77.   } else {
  78.     for (i = 0; i < N; i++) {
  79.       IMAG(C, ldc * i + i) = 0.0;
  80.     }
  81.   }
  82.  
  83.   if (alpha_real == 0.0 && alpha_imag == 0.0)
  84.     return;
  85.  
  86.   if (uplo == CblasUpper && trans == CblasNoTrans) {
  87.  
  88.     for (i = 0; i < N; i++) {
  89.  
  90.       /* Cii += alpha Aik conj(Bik) + conj(alpha) Bik conj(Aik) */
  91.       {
  92.     BASE temp_real = 0.0;
  93.     /* BASE temp_imag = 0.0; */
  94.     for (k = 0; k < K; k++) {
  95.       const BASE Aik_real = CONST_REAL(A, i * lda + k);
  96.       const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  97.       /* temp1 = alpha * Aik */
  98.       const BASE temp1_real = alpha_real * Aik_real - alpha_imag * Aik_imag;
  99.       const BASE temp1_imag = alpha_real * Aik_imag + alpha_imag * Aik_real;
  100.       const BASE Bik_real = CONST_REAL(B, i * ldb + k);
  101.       const BASE Bik_imag = CONST_IMAG(B, i * ldb + k);
  102.       temp_real += temp1_real * Bik_real + temp1_imag * Bik_imag;
  103.     }
  104.  
  105.     REAL(C, i * ldc + i) += 2 * temp_real;
  106.     IMAG(C, i * ldc + i) = 0.0;
  107.       }
  108.  
  109.       /* Cij += alpha Aik conj(Bjk) + conj(alpha) Bik conj(Ajk) */
  110.       for (j = i + 1; j < N; j++) {
  111.     BASE temp_real = 0.0;
  112.     BASE temp_imag = 0.0;
  113.     for (k = 0; k < K; k++) {
  114.       const BASE Aik_real = CONST_REAL(A, i * lda + k);
  115.       const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  116.       /* temp1 = alpha * Aik */
  117.       const BASE temp1_real = alpha_real * Aik_real - alpha_imag * Aik_imag;
  118.       const BASE temp1_imag = alpha_real * Aik_imag + alpha_imag * Aik_real;
  119.       const BASE Bik_real = CONST_REAL(B, i * ldb + k);
  120.       const BASE Bik_imag = CONST_IMAG(B, i * ldb + k);
  121.  
  122.       const BASE Ajk_real = CONST_REAL(A, j * lda + k);
  123.       const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
  124.       /* temp2 = alpha * Ajk */
  125.       const BASE temp2_real = alpha_real * Ajk_real - alpha_imag * Ajk_imag;
  126.       const BASE temp2_imag = alpha_real * Ajk_imag + alpha_imag * Ajk_real;
  127.       const BASE Bjk_real = CONST_REAL(B, j * ldb + k);
  128.       const BASE Bjk_imag = CONST_IMAG(B, j * ldb + k);
  129.  
  130.       /* Cij += alpha * Aik * conj(Bjk) + conj(alpha) * Bik * conj(Ajk) */
  131.       temp_real += ((temp1_real * Bjk_real + temp1_imag * Bjk_imag)
  132.             + (Bik_real * temp2_real + Bik_imag * temp2_imag));
  133.       temp_imag += ((temp1_real * (-Bjk_imag) + temp1_imag * Bjk_real)
  134.             + (Bik_real * (-temp2_imag) + Bik_imag * temp2_real));
  135.     }
  136.     REAL(C, i * ldc + j) += temp_real;
  137.     IMAG(C, i * ldc + j) += temp_imag;
  138.       }
  139.     }
  140.  
  141.   } else if (uplo == CblasUpper && trans == CblasConjTrans) {
  142.  
  143.     for (k = 0; k < K; k++) {
  144.       for (i = 0; i < N; i++) {
  145.     BASE Aki_real = CONST_REAL(A, k * lda + i);
  146.     BASE Aki_imag = CONST_IMAG(A, k * lda + i);
  147.     BASE Bki_real = CONST_REAL(B, k * ldb + i);
  148.     BASE Bki_imag = CONST_IMAG(B, k * ldb + i);
  149.     /* temp1 = alpha * conj(Aki) */
  150.     BASE temp1_real = alpha_real * Aki_real - alpha_imag * (-Aki_imag);
  151.     BASE temp1_imag = alpha_real * (-Aki_imag) + alpha_imag * Aki_real;
  152.     /* temp2 = conj(alpha) * conj(Bki) */
  153.     BASE temp2_real = alpha_real * Bki_real - alpha_imag * Bki_imag;
  154.     BASE temp2_imag = -(alpha_real * Bki_imag + alpha_imag * Bki_real);
  155.  
  156.     /* Cii += alpha * conj(Aki) * Bki + conj(alpha) * conj(Bki) * Aki */
  157.     {
  158.       REAL(C, i * lda + i) += 2 * (temp1_real * Bki_real - temp1_imag * Bki_imag);
  159.       IMAG(C, i * lda + i) = 0.0;
  160.     }
  161.  
  162.     for (j = i + 1; j < N; j++) {
  163.       BASE Akj_real = CONST_REAL(A, k * lda + j);
  164.       BASE Akj_imag = CONST_IMAG(A, k * lda + j);
  165.       BASE Bkj_real = CONST_REAL(B, k * ldb + j);
  166.       BASE Bkj_imag = CONST_IMAG(B, k * ldb + j);
  167.       /* Cij += alpha * conj(Aki) * Bkj + conj(alpha) * conj(Bki) * Akj */
  168.       REAL(C, i * lda + j) += (temp1_real * Bkj_real - temp1_imag * Bkj_imag)
  169.           + (temp2_real * Akj_real - temp2_imag * Akj_imag);
  170.       IMAG(C, i * lda + j) += (temp1_real * Bkj_imag + temp1_imag * Bkj_real)
  171.           + (temp2_real * Akj_imag + temp2_imag * Akj_real);
  172.     }
  173.       }
  174.     }
  175.  
  176.   } else if (uplo == CblasLower && trans == CblasNoTrans) {
  177.  
  178.     for (i = 0; i < N; i++) {
  179.  
  180.       /* Cij += alpha Aik conj(Bjk) + conj(alpha) Bik conj(Ajk) */
  181.  
  182.       for (j = 0; j < i; j++) {
  183.     BASE temp_real = 0.0;
  184.     BASE temp_imag = 0.0;
  185.     for (k = 0; k < K; k++) {
  186.       const BASE Aik_real = CONST_REAL(A, i * lda + k);
  187.       const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  188.       /* temp1 = alpha * Aik */
  189.       const BASE temp1_real = alpha_real * Aik_real - alpha_imag * Aik_imag;
  190.       const BASE temp1_imag = alpha_real * Aik_imag + alpha_imag * Aik_real;
  191.       const BASE Bik_real = CONST_REAL(B, i * ldb + k);
  192.       const BASE Bik_imag = CONST_IMAG(B, i * ldb + k);
  193.  
  194.       const BASE Ajk_real = CONST_REAL(A, j * lda + k);
  195.       const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
  196.       /* temp2 = alpha * Ajk */
  197.       const BASE temp2_real = alpha_real * Ajk_real - alpha_imag * Ajk_imag;
  198.       const BASE temp2_imag = alpha_real * Ajk_imag + alpha_imag * Ajk_real;
  199.       const BASE Bjk_real = CONST_REAL(B, j * ldb + k);
  200.       const BASE Bjk_imag = CONST_IMAG(B, j * ldb + k);
  201.  
  202.       /* Cij += alpha * Aik * conj(Bjk) + conj(alpha) * Bik * conj(Ajk) */
  203.       temp_real += ((temp1_real * Bjk_real + temp1_imag * Bjk_imag)
  204.             + (Bik_real * temp2_real + Bik_imag * temp2_imag));
  205.       temp_imag += ((temp1_real * (-Bjk_imag) + temp1_imag * Bjk_real)
  206.             + (Bik_real * (-temp2_imag) + Bik_imag * temp2_real));
  207.     }
  208.     REAL(C, i * ldc + j) += temp_real;
  209.     IMAG(C, i * ldc + j) += temp_imag;
  210.       }
  211.  
  212.       /* Cii += alpha Aik conj(Bik) + conj(alpha) Bik conj(Aik) */
  213.       {
  214.     BASE temp_real = 0.0;
  215.     /* BASE temp_imag = 0.0; */
  216.     for (k = 0; k < K; k++) {
  217.       const BASE Aik_real = CONST_REAL(A, i * lda + k);
  218.       const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  219.       /* temp1 = alpha * Aik */
  220.       const BASE temp1_real = alpha_real * Aik_real - alpha_imag * Aik_imag;
  221.       const BASE temp1_imag = alpha_real * Aik_imag + alpha_imag * Aik_real;
  222.       const BASE Bik_real = CONST_REAL(B, i * ldb + k);
  223.       const BASE Bik_imag = CONST_IMAG(B, i * ldb + k);
  224.       temp_real += temp1_real * Bik_real + temp1_imag * Bik_imag;
  225.     }
  226.  
  227.     REAL(C, i * ldc + i) += 2 * temp_real;
  228.     IMAG(C, i * ldc + i) = 0.0;
  229.       }
  230.     }
  231.  
  232.   } else if (uplo == CblasLower && trans == CblasConjTrans) {
  233.  
  234.     for (k = 0; k < K; k++) {
  235.       for (i = 0; i < N; i++) {
  236.     BASE Aki_real = CONST_REAL(A, k * lda + i);
  237.     BASE Aki_imag = CONST_IMAG(A, k * lda + i);
  238.     BASE Bki_real = CONST_REAL(B, k * ldb + i);
  239.     BASE Bki_imag = CONST_IMAG(B, k * ldb + i);
  240.     /* temp1 = alpha * conj(Aki) */
  241.     BASE temp1_real = alpha_real * Aki_real - alpha_imag * (-Aki_imag);
  242.     BASE temp1_imag = alpha_real * (-Aki_imag) + alpha_imag * Aki_real;
  243.     /* temp2 = conj(alpha) * conj(Bki) */
  244.     BASE temp2_real = alpha_real * Bki_real - alpha_imag * Bki_imag;
  245.     BASE temp2_imag = -(alpha_real * Bki_imag + alpha_imag * Bki_real);
  246.  
  247.     for (j = 0; j < i; j++) {
  248.       BASE Akj_real = CONST_REAL(A, k * lda + j);
  249.       BASE Akj_imag = CONST_IMAG(A, k * lda + j);
  250.       BASE Bkj_real = CONST_REAL(B, k * ldb + j);
  251.       BASE Bkj_imag = CONST_IMAG(B, k * ldb + j);
  252.       /* Cij += alpha * conj(Aki) * Bkj + conj(alpha) * conj(Bki) * Akj */
  253.       REAL(C, i * lda + j) += (temp1_real * Bkj_real - temp1_imag * Bkj_imag)
  254.           + (temp2_real * Akj_real - temp2_imag * Akj_imag);
  255.       IMAG(C, i * lda + j) += (temp1_real * Bkj_imag + temp1_imag * Bkj_real)
  256.           + (temp2_real * Akj_imag + temp2_imag * Akj_real);
  257.     }
  258.  
  259.     /* Cii += alpha * conj(Aki) * Bki + conj(alpha) * conj(Bki) * Aki */
  260.     {
  261.       REAL(C, i * lda + i) += 2 * (temp1_real * Bki_real - temp1_imag * Bki_imag);
  262.       IMAG(C, i * lda + i) = 0.0;
  263.     }
  264.       }
  265.     }
  266.   } else {
  267.     BLAS_ERROR("unrecognized operation");
  268.   }
  269. }
  270.